Turbulence in non-integer dimensions by fractal Fourier decimation 
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Fractal decimation reduces the effective dimensionality of a flow by keeping only a (randomly 
chosen) set of Fourier modes whose number in a ball of radius k is proportional to k D for large k. 
At the critical dimension D = 4/3 there is an equilibrium Gibbs state with a fc -5 ^ 3 spectrum, as in 
[V. L'vov et al, Phys. Rev. Lett. 89, 064501 (2002)]. Spectral simulations of fractally decimated 
two-dimensional turbulence show that the inverse cascade persists below D — 2 with a rapidly rising 
Kolmogorov constant, likely to diverge as (D — 4/3)" 
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PACS numbers: 47.27.Gs, 05.20.Jj 

In theoretical physics a number of interesting results 
have been obtained by extending the dimension d of space 
from directly relevant values such as 1, 2, 3 to non-integer 
values. Dimensional regularization in field theory jlj and 
the 4 — e expansion in critical phenomena are well- 
known instances. In such approaches, one usually ex- 
pands the solution in terms of Feynman diagrams, each 
of which can be analytically continued to real or complex 
values of d. The same kind of extension can be carried 
out for homogeneous isotropic turbulence but a severe 
difficulty appears then for d < 2: the energy spectrum 
E(k) can become negative in some band of wavenumbers 
k, so that this kind of extension lacks probabilistic real- 
lizability 0. Nevertheless, in Ref. henceforth cited as 
LPP, it is argued that, should there exist an alternative 
realizable way of doing the extension below dimension 
two in which the nonlinearity conserves energy and en- 
strophy, then an interesting phenomenon - to which we 
shall come back - should happen in dimension 4/3. 

For diffusion and phase transitions there is a very 
different way of switching to non-integer dimensions, 
namely to reformulate the problem on a fractal of di- 
mension D (here a capital D will always be a fractal 
dimension) [5|. Are we able to do this for hydrodynam- 
ics? Implementing mass and momentum conservation on 
a fractal is quite a challenged] We discovered a new way 
of fractal decimation in Fourier space, appropriate for 
hydrodynamics. Since, here, we are primarily interested 
in dimensions less than two, we shall do our decimation 
starting from the standard d = 2 case. 

The forced incompressible Navier-Stokes equations for 
the velocity field can be written in abstract notation as 

d t u = B{u,u) + f + Au , (1) 
B(u,u) = -u-Vu + Vp , N = vV 2 , (2) 

where u stands for the velocity field u(xi, x%, t), f for 
f(xi,X2,t), p is the pressure and v the viscosity. The 



velocity u is taken in the space of divergence- less velocity 
fields which are 2-7T periodic in x% and X2, such that u(t = 
0) = Uq. Now, we define a Fourier decimation operator 
Pd on this space of velocity fields: 



If u = 



kez 2 



°Uk : then Pdu = 



kez 2 



TkUk 



Here 0k are random numbers such that 

I 1 with probability hk 
with probability 1 — hk , f 



\k\ 



(3) 



(4) 



To obtain D-dimcnsional dynamics we choose 

h k = C{k/k ) D -' 2 , < D < 2, < C < 1 , (5) 

where k$ is a reference wavenumber; here C = kq = 1. 
All the 0k are chosen independently, except that 8k = 
0-k to preserve Hermitian symmetry. Our fractal deci- 
mation procedure removes at random — but in a time- 
frozen (quenched) way — many modes from the k lattice, 
leaving on average N(k) oc kP active modes in a disk of 
radius k. The randomness in the choice of the decimation 
will be called the disorder^ 

Observe that Pd is a projector, that it commutes 
with the viscous diffusion operator A and that it is self- 
adjoint for the energy (L 2 ) norm, defined as usual as 
||it|| 2 = (l/(27r) 2 ) J \u(x)\ 2 d 2 x, where the integral is 
over a 27r x 27T periodicity square. The conservation of 
energy (by the nonlinear term) for sufficiently smooth so- 
lutions of the Navier-Stokes equation can be expressed 
as (u, B(u,u)) = where (it, w) = (l/(27r) 2 ) J u(x) ■ 
w(x) d 2 x is the L 2 scalar product. 

The decimated Navier-Stokes equation, written for an 
incompressible field v, takes the following form: 



d t v = P D B{v,v) + P D f + P D Av 



(6) 



1 Lattice Boltzmann models may be amenable to fractal decima- 
tion [f|. 



2 A more drastic and non-random decimation is the reduced wave- 
vector set approximation (REWA; see Ref. and references 
therein), in which the number of active mode grows as lnfc, so 
that from our point of view it has dimension D = 0. 



2 



The initial condition is vq = v(t = 0) = PdUq. Thus 
at any later time Pdv = v. Energy is again conserved; 
indeed (v, PrjB(v, v)) = 0, as is seen by moving the self- 
adjoint operator P D to the left hand side of the scalar 
product and using PdV = v. Enstrophy conservation 
for the decimated problem is proved in a smilar way by 
working with the vorticity. 

If, in addition to decimation, we apply a Galerkin trun- 
cation which kills all the modes having wavenumbers be- 
yond a threshold Kq, the surviving modes constitute a 
dynamical system having a finite number of degrees of 
freedom. Such truncated systems with no forcing and no 
viscosity have been studied by Lee, Kraichnan and others 
[|| . Using suitable variables related to the real and imag- 
inary parts of the active modes, the dynamical equations 
may be written as y a = J^p^ A a p 7 ypy 7 . 

For the purely Galerkin-truncated (not decimated) 
case it is well known that the above dynamical sys- 
tem satisfies a Liouville theorem J2 a dy a /dy a = and 
thus conserves volume in phase space. This in turn im- 
plies the existence of (statistically) invariant Gibbs states 
for which the probability is a Gaussian, proportional 
to e -( aE +0 n ) where E — ^2 k \uk\ 2 is the energy and 
Q = "^2 k k 2 \uk\ 2 is the enstrophy. Such Gibbs states, 
called by Kraichnan absolute equilibria, play an impor- 
tant role in his theory of the two-dimensional (2D) inverse 
energy cascade 0. If we now combine inviscid, unforced 
Galerkin truncation and decimation, it is easily checked 
that the Liouville theorem still holds, provided the dec- 
imation preserves Hermitean symmetry. For such Gibbs 
states, and any active mode (9k — 1), one easily checks 
that the mean square energy (|«fc| 2 ) = C'/(a + 13k 2 ), 
where C > does not depend on k. The correspond- 
ing energy spectrum is the mean energy E(k) of modes 
having a wavenumber between k and k + 1. Up to fluc- 
tuations of the disorder, the number of active modes in 
such a shell is Oik ' 1 ). Thus, 



E(k) = 



,D-\ 



a + /3k 2 '' 



(3 > 0, a > -f3, 



(7) 



where various positive constants have been absorbed into 
a new definition of a and (3. An instance is enstrophy 
equipartition: a = (all the modes have the same en- 
strophy), for which the energy spectrum is E(k) oc k D ~ 3 . 
Following LPP, if we now set D = 4/3, we obtain a k~ 5 ^ 3 
spectrum, the spectrum predicted by Kolmogorov's 1941 
scaling theory and extended by Kraichnan to the inverse 
energy cascade of 2D turbulence. Note that such Gibbs 
states are only conditionally Gaussian, for a given dis- 
order. Otherwise, they are highly intermittent, since a 
given high-fc mode will be active only in a small fraction 
of the disorder realizations. We also note that similar 
phenomena have been observed in shell models fioj . 

The form (J7J of the D-dimensional absolute equilib- 
ria also allows for the kind of Bose condensation in the 
gravest modes (here, those with unit wavenumber) found 



by Kraichnan for 2D turbulence. For this the "inverse 
temperature" a must be taken negative, close to its 
minimum realizable value —(3. The arguments used by 
Kraichnan to predict an inverse Kolmogorov fc -5 / 3 en- 
ergy cascade for high-Reynolds number 2D turbulence 
with forcing near an intermediate wavenumber fci n j carry 
over to the decimated case with D < 2. In particular the 
conservation of enstrophy blocks energy transfer to high 
wavenumbers. This in itself does not imply that the en- 
ergy will cascade to wavenumbers smaller than fc; n j , pro- 
ducing a fc-independent energy flux: it might also linger 
around and accumulate near fc; n j . 

It is now our purpose to show that for 4/3 < D < 2, 
when the energy spectrum is prescribed to be E(k) = 
£-5/3 over t ne inertial range, there is a negative energy 
flux T\e, vanishing linearly with D — 4/3 near the criti- 
cal dimension D — 4/3. For this we shall assume that a 
key feature of the two-dimensional energy cascade carries 
over to lower dimensions, namely the existence of scaling 
solutions with local (in Fourier space) dynamics, so that 
the energy transfer is dominated by triads of wavenum- 
bers with comparable magnitudes. Let us now decom- 
pose the energy inertial range into bands of fixed relative 
width, say one octave, delimited by the wavenumbers 
2°, 2 1 , 2 2 , etc. Because of locality there is much intra- 
band dynamics but, of course, interband interactions are 
needed to obtain an energy flux. Pure intraband dy- 
namics (with no forcing and dissipation) would lead to 
thermalization. For dimensional reasons, thermalization 
and interband transfer have the same time scale, namely 
the eddy turnover time k~ 3 ^ 2 E~ 1 ^ 2 (k). 

To get a handle on the combined intraband and inter- 
band dynamics we perform a thermodynamic thought ex- 
periment in which we artificially separate them in time. 
In the first phase, starting from a k~ 5 ^ 3 spectrum we 
prevent the various bands from interacting by introduc- 
ing (impenetrable) interband barriers at their edges. In 
each band, the modes will then thermalize and achieve a 
Gibbs state with a spectrum (0 in which a and (3 are de- 
termined by the constraints that the total band energy 
and enstrophy be the same as for the —5/3 spectrum. 
For example, in the first band this gives the constraints 
(n = for the energy and n = 2 for the enstrophy) 



dkk n 



k D - x l (a + f3k 2 ) - k 



-5/3 



0. 



(8) 



In general this constitutes a system of transcendental 
equations for the parameters a and f3 that can only be 
solved numerically. This is illustrated in Fig. [T] for the 
2D case. We observe that in 2D the absolute equilibrium 
spectrum and the —5/3 spectrum are very close to each 
other. Specifically, in 2D the absolute equilibrium spec- 
trum exceeds the —5/3 spectrum by about 10% at any 
lower band edge and by about 5% at any upper band 
edge. Of course, as we approach the critical dimension 
D = 4/3 the discrepancy goes to zero and can easily be 
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FIG. 1. (Color online) The k~ 5 ^ 3 spectrum (continuous) and 
the associated 2D absolute equilibrium with the same energy 
and enstrophy in the first octave (dashed). 



calculated pertubatively in D — 4/3. In the second phase 
of our thought experiment, we remove just one of the bar- 
riers between two adjacent bands, say, the barrier at 2 1 . 
A new thermalization leads then to an absolute equilib- 
rium in the band [2°, 2 2 ], which again, can be easily cal- 
culated. In 2D, before the removal, the energy between 
2° and 2 1 was 0.555. After the new thermalization, this 
energy is found to have increased by 0.00551. Thus en- 
ergy has been transferred from the upper band [2 1 , 2 2 ] to 
the lower band [2°, 2 1 ]. Close to D — 4/3, we can again 
apply elementary perturbation techniques and obtain for 
the upper-to-lower-band energy transfer 0.009(.D — 4/3) 
to leading order. Our thermodynamic thought experi- 
ment thus suggests that the energy flux vanishes linearly 
with D — 4/3, being negative above the critical dimen- 
sion, which implies an inverse cascade. In the K41 inertial 
range, the energy spectrum and the energy flux He are 
related by E{k) = CKoi|n B | 2 / 3 fc~ 5 / 3 , where C Ko i is the 
Kolmogorov constant, we infer that the Kolmogorov con- 
stant diverges as (D — 4/3)~ 2 / 3 . A closure calculation of 
eddy-damped quasi-normal Markovian (EDQNM) type 
also predicts a divergence with a —2/3 exponent. 

Kraichnan's ideas about the inverse cascade in 2D got 
growing support a few years later from direct numerical 
simulations, which eventually achieved the resolution of 
32, 768 2 modes " 
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As to our idea about the robust- 
ness of the inverse cascade and the growth of the Kol- 
mogorov constant when lowering the dimension D, some 
support can be already provided, using a D-dimensional 
decimated variant of spectral direct numerical simula- 
tion: First one generates an instance of the disorder, that 
is the list of active and inactive Fourier modes; then, one 
applies standard time marching algorithms and, at each 
time step, sets to zero all inactive modes. In addition to 
the well-known difficulties of simulating 2D turbulence 
(see e.g. 11| and references therein), there are new diffi- 



culties. 

A few words about the numerical implementation. We 



integrate the decimated Navier-Stokes equation ([6]) in 
vorticity representation. Instead of using as damping the 
viscous operator A = vA (where A = V 2 is the Lapla- 
cian), we use 



A 



-is A 



+2 



//A" 



v > 0, \i > 0, 



(9) 



whose Fourier symbol is —vk A — pik~ 4 . In other words, 
we use hyperviscosity to avoid wasting resolution on the 
enstrophy cascade and large-scale friction to prevent an 
accumulation of energy on the gravest modes and thus 
allow eventual convergence to a statistical steady state. 
The results reported here have a resolution of N = 3072 
collocation points in the two coordinates. Time march- 
ing is done by an Adams-Bashforth scheme combined 
with exponential time difference (ETD) [lij with a time 
step between 5 x 10~ 5 and 10~ 4 , depending on dimen- 
sion. Energy injection at the rate e is done in a band 
of width three around k m j — 319 by adding to the time- 
rate-of change of the Fourier amplitude of the vorticity a 
term proportional to the inverse of its complex conjugate 
13]. This allows a fc-independent and time-independent 
energy injection. As D is decreased the amplitude of this 
forcing is increased to keep the total energy injection on 
active modes fixed at e = 0.01. The damping parameters 
are v = 10~ n and /i = 0.1. Runs are done concurrently 
for different values of D on a high-performance cluster 
at the Weizmann Institute and take typically a few thou- 
sand hours of CPU per run to achieve a statistical steady 
state. 

Energy spectra are obtained by angular averages over 
Fourier-space shells of unit width 
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(10) 



K<k<K+l 



where the v(k) are the Fourier coefficients of the solution 
of the decimated Navier-Stokes equation ([6]). We also 
need the energy flux He(K) through wavenumber K due 
to nonlinear transfer, defined as 



U E (K) = «*(*:) -NL(fe), 



(11) 



k<K 



where NL(fe) denotes the set of Fourier coefficients of 
the nonlinear term Pe>B(v,v) in the decimated Navier- 
Stokes equation (|6]) and the asterisk denotes complex 
conjugation. E{k) and He(K) are mostly insensitive to 
the disorder realization. 

Figs. [2] and [3] (inset) show the steady-state compen- 
sated energy spectra k 5 ^ 3 E(k) and the energy fluxes 
Hsik), for various values of D between 2 and 1.5, re- 
spectively. Both are quite flat, over a significant range 
of k values, evidence that D-dimensional forced turbu- 
lence, Fourier decimated down from the two-dimensional 
case, preserves the key feature of two-dimensional turbu- 
lence of having an inverse cascade that follows the —5/3 
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FIG. 2. (Color online) Compensated steady-state spectra for 
D = 2.0, 1.9, 1.8, 1.7, 1.6, 1.5 from bottom to top with spikes 
at injection. The inset shows the dependence on D of the 
plateau of the compensated spectra, as an average over the 
interval between vertical dashed lines (with standard devia- 
tion error bars). 
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FIG. 3. (Color online) Dependence of the Kolmogorov con- 
stant on D. The lowest value, at D = 2, is about 5. The inset 
shows the energy flux normalized by the energy injection s 
for the same values of D as in Fig. [2] 



law. Note that the inertial range (the flat region of the 
compensated energy spectrum) shrinks as the dimension 
D decreases. The absolute value of the energy flux is 
about 80% of the energy injection e for D = 2, but drops 
to less than 50% for D — 1.5. Indeed, as the dimen- 
sion D is lowered, there are fewer and fewer pairs of ac- 
tive modes in the forcing band, capable through their 
beating interaction of draining the energy into the in- 
frared direction; thus the energy injection will be more 
and more balanced by direct dissipation near injection. 
Preventing this would require a substantial lowering of 
the dissipation which in turn requires a substantial in- 



crease in the resolution at the high-fc end. Anyway, the 
fact that the flux III^I becomes substantially lower than 
injection does not prevent us from calculating the Kol- 
mogorov constant, given (in terms of plateau values) by 
CkoI = fc 5/3 £(fc)/ (\U E (k)\ 2 ^). Figured shows the vari- 
ation of the Kolmogorov constant with dimension. When 
lowering the dimension from 2 to 1.5, a combined effect 
of a rise in the compensated spectrum and a drop in flux 
yields a monotonic growth of about a factor ten in the 
Kolmogorov constant and a substantial growth of errors 
due to fluctuations within the averaging interval. Prob- 
ing the conjectured divergence by moving closer to the 
critical point D = 4/3 would require much higher resolu- 
tion. A state-of-the-art 16, 384 2 simulation of sufficient 
length might shed light. 

We finally observe that the fractal Fourier decimation 
procedure can be started from any integer dimension and 
can be applied to a large class of problems in hydrody- 
namics and beyond. It could be interesting, for example, 
to investigate how it affects the dissipative anomaly of 
shock-dominated compressible flow. 
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